In [9]:
import string
import copy
import scipy
import Tkinter, tkFileDialog
import numpy as np
import pandas as pd
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
import os
import sys
import datetime
import matplotlib.units
import re
from numba import jit,int32
import time
sys.path.append(os.path.abspath("C:\Users\Scherer Lab E\Documents\GitHub\Python_Data_Analysis"))
import common_functions
%matplotlib inline

Other Data Set

In [10]:
from IPython.display import Image
os.chdir("C:\Users\Scherer Lab E\Documents\IPython Notebooks\SPIFF Manuscript\Data")
Image(filename="06191317_example_image.png")
Out[10]:
In [11]:
os.chdir("C:\Users\Scherer Lab E\Documents\IPython Notebooks\SPIFF Manuscript\Data")
Image(filename="06191317_particle_closeup.png")
Out[11]:

Trackpy

Use Trackpy with different window sizes to get different SPIFFs of this data. Chose windows of 7, 9, 11, 13, and 15 pixels

In [2]:
import pims
import trackpy as tp
from astropy.io import fits
In [3]:
frames = pims.ImageSequence("E:\Pat's Projects\Constant Dragging\Exp06191301\Mov_06191306 - 200fps 350x2560\*.tif")
In [57]:
'''Localize the image data using several windows'''
localized_data = []
for i in range(5):
    f_other_data = tp.batch(frames, 7+i*2, separation=5)
    localized_data.append([f_other_data, 7+i*2])

Frame 974: 4061 features
In [28]:
'''Checkpoint after localization'''
import cPickle
pik = open('Ana_16041501_localized_positions.pkl', 'w')
cPickle.dump(localized_data, pik)
pik.close()
In [29]:
import cPickle
pik = open('Ana_16041501_localized_positions.pkl', 'r')
localized_data = cPickle.load(pik)
pik.close()
In [30]:
'''Link the image data for each window'''
tracked_data = []
for frame, window in localized_data:
    tracked = tp.link_df(frame, 10, memory=1)
    tracked_data.append([tracked, window])

Frame 974: 4061 trajectories present
In [31]:
'''Checkpoint after linking'''
import cPickle
pik = open('Ana_16041501_linked_positions.pkl', 'w')
cPickle.dump(tracked_data, pik)
pik.close()
In [32]:
'''Load the linked data'''
import cPickle
pik = open('Ana_16041501_linked_positions.pkl', 'r')
tracked_data = cPickle.load(pik)
pik.close()
In [37]:
def disp_calc_en(grp, tau):
    if len(grp) < tau:
        return pd.Series([np.nan])
    #print grp[['x pos', 'y pos']],grp.shift(tau, axis=0)[['x pos', 'y pos']]
    disp = ((grp[['x pos', 'y pos']] - grp.shift(tau, axis=0)[['x pos', 'y pos']])**2).sum(skipna=False,axis=1)
    #if grp.name > 5000:
        #print grp.name
    return disp

def calc_en_msd(df, percent=0.3):
    en_msd = []
    pivot = pd.pivot_table(df, values=['x pos','y pos'], index='frame', columns='track id')
    for tau in range(int(max(df['frame'])*percent)):
        del_xy = (pivot - pivot.shift(tau+1, axis=0))**2
        en_msd.append((del_xy['x pos'] + del_xy['y pos']).stack().mean())
    return en_msd
In [61]:
#%pdb on
temp_track = tracked_data[0][0]
temp_track = temp_track.rename(columns={'x': 'x pos', 'y': 'y pos', 'particle': 'track id' })
temp_track['track id'] = temp_track['track id'].astype('int64')
temp_track_sorted = temp_track.sort_values(['frame', 'track id'])
#temp_track.groupby('track id').apply(lambda x: x.name)

Calculated the ensemble averaged MSD on the first 3rd of the data.

In [38]:
'''Calculate the ensemble MSD for each window (for only first third of tau)'''
en_msd = []
for tracks, window in tracked_data:
    temp_track = tracks.rename(columns={'x': 'x pos', 'y': 'y pos', 'particle': 'track id' })
    en_msd.append([calc_en_msd(temp_track),window])
    print 'done with window '+str(window)

done with window 7
done with window 9
done with window 11
done with window 13
done with window 15
In [43]:
'''Checkpoint after calculating ensemble MSD'''
import cPickle
pik = open('Ana_16041501_ensemble_msd.pkl', 'w')
cPickle.dump(en_msd, pik)
pik.close()
In [44]:
'''Load the Checkpointed MSD Data'''
import cPickle
os.chdir("C:\Users\Scherer Lab E\Documents\IPython Notebooks\SPIFF Manuscript\Data")
pik = open('Ana_16041501_ensemble_msd.pkl', 'r')
en_msd = cPickle.load(pik)
pik.close()
In [45]:
for num, entry in enumerate(tracked_data):
    df, window = entry
    plt.figure(figsize=[10,4])
    plt.subplot(121)
    plt.hist2d(df.x%1, df.y%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    
    len_msd = len(en_msd[num][0])+1
    ax2 = plt.subplot(122)
    plt.plot(range(1,len_msd), en_msd[num][0], '-ob')
    pwr_law = lambda x, param: param*x
    fit_params = scipy.optimize.curve_fit(pwr_law, np.arange(1,len_msd)[-200:], en_msd[num][0][-200:])
    plt.plot(np.arange(1,len_msd), pwr_law(np.arange(1,len_msd), fit_params[0][0]), 'g')
    plt.loglog()
    ax2.text(0.5, 0.2, str(fit_params[0][0])+'t', transform=ax2.transAxes)
    plt.title('MSD and Power Law for window='+str(window))
    plt.ylabel('MSD (Pixels$^2$)')
    plt.xlabel('Lag Time (Frames)')
    plt.tight_layout()
    plt.show()

Find Residuals of the MSD Fit

To see if the SPIFF affects the MSD at all we can subtract the power law fit from the MSD to see if the change in the MSD with the SPIFF is more sublte.

In [46]:
for num, entry in enumerate(tracked_data):
    df, window = entry
    plt.figure(figsize=[14,4])
    plt.subplot(131)
    plt.hist2d(df.x%1, df.y%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    
    ax2 = plt.subplot(132)
    plt.plot(range(1,len_msd), en_msd[num][0], '-ob')
    pwr_law = lambda x, param: param*x
    fit_params = scipy.optimize.curve_fit(pwr_law, np.arange(1,len_msd)[-240:], en_msd[num][0][-240:])
    plt.plot(np.arange(1,len_msd), pwr_law(np.arange(1,len_msd), fit_params[0][0]), 'g')
    plt.loglog()
    ax2.text(0.5, 0.2, str(round(fit_params[0][0],5))+'t', transform=ax2.transAxes)
    plt.title('MSD and Power Law for window='+str(window))
    plt.ylabel('MSD (Pixels$^2$)')
    plt.xlabel('Lag Time (Frames)')
    
    len_msd = len(en_msd[num][0])+1
    ax3 = plt.subplot(133)
    #plt.plot(range(1,len_msd), en_msd[num][0], '-ob')
    pwr_law = lambda x, param: param*x
    fit_params = scipy.optimize.curve_fit(pwr_law, np.arange(1,len_msd)[-240:], en_msd[num][0][-240:])
    resid = en_msd[num][0] - pwr_law(np.arange(1,len_msd), fit_params[0][0])
    plt.plot(np.arange(1,len_msd), resid, 'b')
    #plt.loglog()
    ax3.text(0.2, 0.2, str(round(fit_params[0][0],5))+'t', transform=ax3.transAxes)
    plt.title('MSD - Power Law Fit for window='+str(window))
    plt.ylabel('MSD (Pixels$^2$)')
    plt.xlabel('Lag Time (Frames)')
    
    plt.show()

SPIFF Correction

In [47]:
def spiff_correction(df):
    spiff_cent = df.copy()
    df_translated = df[['x', 'y']] + 0.5
    spiff_cent[['x', 'y']] = df_translated[['x', 'y']] % 1
    
    # Correct x > 0.5 values
    spiff_x_plus = spiff_cent[spiff_cent.x >= 0.5]
    spiff_x_plus_percentile = spiff_x_plus.x.rank()/len(spiff_x_plus)
    spiff_x_plus_corrected = 0.5 + spiff_x_plus_percentile * 0.5
    
    spiff_x_minus = spiff_cent[spiff_cent.x < 0.5]
    spiff_x_minus_percentile = spiff_x_minus.x.rank(ascending=False)/len(spiff_x_minus)
    spiff_x_minus_corrected = 0.5 - spiff_x_minus_percentile * 0.5
    
    spiff_y_plus = spiff_cent[spiff_cent.y >= 0.5]
    spiff_y_plus_percentile = spiff_y_plus.y.rank()/len(spiff_y_plus)
    spiff_y_plus_corrected = 0.5 + spiff_y_plus_percentile * 0.5
    
    spiff_y_minus = spiff_cent[spiff_cent.y < 0.5]
    spiff_y_minus_percentile = spiff_y_minus.y.rank(ascending=False)/len(spiff_y_minus)
    spiff_y_minus_corrected = 0.5 - spiff_y_minus_percentile * 0.5
    
    x_corrected = pd.concat([spiff_x_plus_corrected, spiff_x_minus_corrected], axis=0)
    y_corrected = pd.concat([spiff_y_plus_corrected, spiff_y_minus_corrected], axis=0)
    
    xy_corrected = pd.DataFrame()#pd.DataFrame([x_corrected, y_corrected], columns=['x corr', 'y corr'])
    xy_corrected['x corr'] = x_corrected.sort_index()
    xy_corrected['y corr'] = y_corrected.sort_index()

    xy_corrected['x corr'] += np.floor(df_translated['x'])
    xy_corrected['y corr'] += np.floor(df_translated['y'])
    xy_corrected = xy_corrected - 0.5
    
    spiff_cent[['x', 'y']] = xy_corrected[['x corr', 'y corr']]
    
    return spiff_cent
In [48]:
corr_tracked_data = []
for tracks, window in tracked_data:
    tracks_temp = spiff_correction(tracks)
    corr_tracked_data.append([tracks_temp, window])
In [49]:
for num, [tracks, window] in enumerate(tracked_data):
    tracks_temp = spiff_correction(tracks)
    
    plt.figure(figsize=[10,4])
    plt.subplot(121)
    plt.hist2d((tracks.x+0.5)%1, (tracks.y+0.5)%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    
    plt.subplot(122)
    plt.hist2d((tracks_temp.x+0.5)%1, (tracks_temp.y+0.5)%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('Corrected SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    plt.show()

MSD of SPIFF Corrected Data

In [50]:
'''Calculate the ensemble MSD for each window (for only first third of tau)'''
corr_en_msd = []
for tracks, window in corr_tracked_data:
    temp_track = tracks.rename(columns={'x': 'x pos', 'y': 'y pos', 'particle': 'track id' })
    corr_en_msd.append([calc_en_msd(temp_track),window])
    print 'done with window '+str(window)

done with window 7
done with window 9
done with window 11
done with window 13
done with window 15
In [51]:
for num, entry in enumerate(corr_tracked_data):
    df, window = entry
    plt.figure(figsize=[10,4])
    plt.subplot(121)
    plt.hist2d(df.x%1, df.y%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    
    len_msd = len(corr_en_msd[num][0])+1
    ax2 = plt.subplot(122)
    plt.plot(range(1,len_msd), corr_en_msd[num][0], '-ob')
    pwr_law = lambda x, param: param*x
    fit_params = scipy.optimize.curve_fit(pwr_law, np.arange(1,len_msd)[-200:], corr_en_msd[num][0][-200:])
    plt.plot(np.arange(1,len_msd), pwr_law(np.arange(1,len_msd), fit_params[0][0]), 'g')
    plt.loglog()
    ax2.text(0.5, 0.2, str(fit_params[0][0])+'t', transform=ax2.transAxes)
    plt.title('MSD and Power Law for window='+str(window))
    plt.ylabel('MSD (Pixels$^2$)')
    plt.xlabel('Lag Time (Frames)')
    plt.tight_layout()
    plt.show()

SPIFF Correction Revision

Rewrote the SPIFF correction funciton to be a lot faster and much simpler to understand.

In [52]:
def spiff_correction(df):
    df_copy = df.copy()
    spiff = df_copy[['x','y']] % 1
    spiff['spiff_x_corr'] = spiff.x.rank()/len(spiff)
    spiff['spiff_y_corr'] = spiff.y.rank()/len(spiff)
    df_copy['x'] = np.floor(df_copy.x) + spiff.spiff_x_corr
    df_copy['y'] = np.floor(df_copy.y) + spiff.spiff_y_corr
    
    return df_copy
In [53]:
for num, [tracks, window] in enumerate(tracked_data):
    tracks_temp = spiff_correction(tracks)
    
    plt.figure(figsize=[10,4])
    plt.subplot(121)
    plt.hist2d((tracks.x)%1, (tracks.y)%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    
    plt.subplot(122)
    plt.hist2d((tracks_temp.x)%1, (tracks_temp.y)%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('Corrected SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    
    print 'Average particles per frame = '+str(df.groupby('frame').apply(len).mean())
    
    plt.show()

Average particles per frame = 4065.69538462

Average particles per frame = 4065.69538462

Average particles per frame = 4065.69538462

Average particles per frame = 4065.69538462

Average particles per frame = 4065.69538462

Calculate the MSD with the new SPIFF Correction

In [54]:
corr_tracked_data = []
for tracks, window in tracked_data:
    tracks_temp = spiff_correction(tracks)
    corr_tracked_data.append([tracks_temp, window])
In [55]:
'''Calculate the ensemble MSD for each window (for only first third of tau)'''
corr_en_msd = []
for tracks, window in corr_tracked_data:
    temp_track = tracks.rename(columns={'x': 'x pos', 'y': 'y pos', 'particle': 'track id' })
    corr_en_msd.append([calc_en_msd(temp_track),window])
    print 'done with window '+str(window)

done with window 7
done with window 9
done with window 11
done with window 13
done with window 15

The MSDs from this new SPIFF correction algorithm produces the same MSDs as with the old ones. I also tried to do a radially symmetric SPIFF correction but it did not work well at all for the square shaped SPIFFs shown in this notebook.

In [56]:
for num, entry in enumerate(corr_tracked_data):
    df, window = entry
    plt.figure(figsize=[10,4])
    plt.subplot(121)
    plt.hist2d(df.x%1, df.y%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    
    len_msd = len(corr_en_msd[num][0])+1
    ax2 = plt.subplot(122)
    plt.plot(range(1,len_msd), corr_en_msd[num][0], '-ob')
    pwr_law = lambda x, param: param*x
    fit_params = scipy.optimize.curve_fit(pwr_law, np.arange(1,len_msd)[-200:], corr_en_msd[num][0][-200:])
    plt.plot(np.arange(1,len_msd), pwr_law(np.arange(1,len_msd), fit_params[0][0]), 'g')
    plt.loglog()
    ax2.text(0.5, 0.2, str(fit_params[0][0])+'t', transform=ax2.transAxes)
    plt.title('MSD and Power Law for window='+str(window))
    plt.ylabel('MSD (Pixels$^2$)')
    plt.xlabel('Lag Time (Frames)')
    plt.tight_layout()
    
    print 'Average particles per frame = '+str(df.groupby('frame').apply(len).mean())
    
    plt.show()

Average particles per frame = 4213.37128205

Average particles per frame = 4177.21128205

Average particles per frame = 4143.63692308

Average particles per frame = 4109.09641026

Average particles per frame = 4065.69538462

Summary of SPIFF Correction and MSD

In [60]:
for num, entry in enumerate(corr_tracked_data):
    df, window = tracked_data[num]
    plt.figure(figsize=[8,6])
    plt.subplot(221)
    plt.hist2d(df['x']%1, df['y']%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    
    len_msd = len(en_msd[num][0])+1
    ax2 = plt.subplot(222)
    plt.plot(range(1,len_msd), en_msd[num][0], '-ob')
    pwr_law = lambda x, param: param*x
    fit_params = scipy.optimize.curve_fit(pwr_law, np.arange(1,len_msd)[-200:], en_msd[num][0][-200:])
    plt.plot(np.arange(1,len_msd), pwr_law(np.arange(1,len_msd), fit_params[0][0]), 'g')
    plt.loglog()
    ax2.text(0.5, 0.2, str(fit_params[0][0])+'t', transform=ax2.transAxes)
    plt.title('MSD and Power Law for window='+str(window))
    plt.ylabel('MSD (Pixels$^2$)')
    plt.xlabel('Lag Time (Frames)')
    plt.tight_layout()
    
    df, window = entry
    plt.subplot(223)
    plt.hist2d(df['x']%1, df['y']%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('Corrected SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    
    len_msd = len(corr_en_msd[num][0])+1
    ax4 = plt.subplot(224)
    plt.plot(range(1,len_msd), corr_en_msd[num][0], '-ob')
    pwr_law = lambda x, param: param*x
    fit_params = scipy.optimize.curve_fit(pwr_law, np.arange(1,len_msd)[-200:], corr_en_msd[num][0][-200:])
    plt.plot(np.arange(1,len_msd), pwr_law(np.arange(1,len_msd), fit_params[0][0]), 'g')
    plt.loglog()
    ax4.text(0.5, 0.2, str(fit_params[0][0])+'t', transform=ax4.transAxes)
    plt.title('Corrected MSD and Power Law for window='+str(window))
    plt.ylabel('MSD (Pixels$^2$)')
    plt.xlabel('Lag Time (Frames)')
    plt.tight_layout()
    
    print 'Average particles per frame = '+str(df.groupby('frame').apply(len).mean())
    
    plt.show()

Average particles per frame = 4213.37128205

Average particles per frame = 4177.21128205

Average particles per frame = 4143.63692308

Average particles per frame = 4109.09641026

Average particles per frame = 4065.69538462
In [63]:
import cPickle
f = open("Ana_16041501_spiff_results.pkl", 'w')
cPickle.dump([tracked_data, en_msd, corr_tracked_data, corr_en_msd], f)
f.close()
In [66]:
tracked_data_slim = []
corr_tracked_data_slim = []
for num, [df, window] in enumerate(tracked_data):
    tmp = df[['x','y','frame','particle']]
    tracked_data_slim.append([tmp, window])
    
    tmp_2 = corr_tracked_data[num][0][['x','y','frame','particle']]
    corr_tracked_data_slim.append([tmp_2, window])
In [69]:
import cPickle
f = open("Ana_16041501_spiff_results.pkl", 'w')
cPickle.dump([tracked_data_slim, en_msd, corr_tracked_data_slim, corr_en_msd], f)
f.close()
In [68]:
pwd
Out[68]:
u'C:\\Users\\Scherer Lab E\\Documents\\IPython Notebooks\\SPIFF Manuscript\\Data'
In [ ]: